/**
 * @file gauss_elimination.cpp
 * @author shanyoumu (22025028@zju.edu.cn)
 * @brief 高斯消去法的实现
 * 程序的bug：在数据类型为double时，不小心加入了int
 * @version 0.1
 * @date 2021-03-14
 * 
 * @copyright Copyright (c) 2021
 * 
 */

#include <vector>
#include <iostream>
#include <cmath>
#define epsilon 1e-16
/**
 * @brief 输出矩阵从k+1行开始到最后一行的绝对值最大的主元编号
 * 
 * @param A 输入的增广矩阵
 * @param k 开始比较的行数
 * @return int 最大元所在的行号
 */
int max_element(std::vector<std::vector<double>> &A, int k)
{
    int result = k;
    for(int i=k+1; i<A.size();i++)
    {
        if(fabs(A[result][result])<fabs(A[i][i]))
            result = i;
    }
    return result;
}

/**
 * @brief 高斯消去法的主要部分
 * 
 * @param A 输入的增广矩阵
 * @return std::vector<double> 求解方程的结果
 */
std::vector<double> Gauss_elimination(std::vector<std::vector<double>> &A)
{
    double detA = 1;
    std::vector<double> B(A.size(),0);
    for(int k =0; k<A.size()-1;k++)
    {
        int maxelement = max_element(A,k);
        if(fabs(A[maxelement][maxelement])<epsilon)
        {
            detA = 0;
            break;
        }
        if(maxelement != k)//若最大元不是第k+1行，则把最大元所在行与第k+1行交换
        {
            detA =-1 *detA;
            for(int i =k; i<A.size()+1;i++)
            {
                double temp = A[maxelement][i];
                A[maxelement][i] = A[k][i];
                A[k][i] = temp;
            }
        }
        for(int i =k+1;i<A.size();i++)
        {
            A[i][k] = A[i][k]/A[k][k];//用来存储乘积因子
            for(int j = k+1; j<A.size()+1; j++)
            {
                A[i][j] = A[i][j] - A[i][k]*A[k][j];
            }
        }
        detA = A[k][k]*detA;
    }
    if(fabs(A[A.size()-1][A.size()-1])<epsilon)//当最后一行的主元过小时
        return B;
    else{
        detA = A[A.size()-1][A.size()-1] *detA;
        A[A.size()-1][A.size()] = A[A.size()-1][A.size()]/A[A.size()-1][A.size()-1];
        for(int k =A.size()-2;k>=0;k--)
        {
            double tem =0;
            for(int j =k+1;j<A.size();j++)
                tem = tem +A[k][j]*A[j][A.size()];
            A[k][A.size()]= (A[k][A.size()]-tem)/A[k][k];
        }
        for(int k =0;k<A.size();k++)
        {
            B[k] = A[k][A.size()];
        }
        return B;
    }
    

}